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Abstract. To interpret the mean depth of cosmic ray air shower maximum and its disper- 
sion, we parametrize those two observables as functions of the first two moments of the In A 
distribution. We examine the goodness of this simple method through simulations of test 
mass distributions. The application of the parameterization to Pierre Auger Observatory 
data allows one to study the energy dependence of the mean In A and of its variance under 
the assumption of selected hadronic interaction models. We discuss possible implications of 
these dependences in term of interaction models and astrophysical cosmic ray sources. 

Keywords: cosmic ray experiments, ultra high energy cosmic rays 
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1 Introduction 

The most commonly used shower observables for the study of the composition of Ultra High 
Energy Cosmic Rays (UHECR) are the mean value of the depth of shower maximum, (A max ), 
and its dispersion, cr(X max ). Inferring the mass composition from these measurements is 
subject to some level of uncertainty. This is because their conversion to mass relies on the 
use of shower simulation codes which include the assumption of a hadronic interaction model. 
The various interaction models [1] have in common the ability to fit lower energy accelerator 
data. However, different physical assumptions are used to extrapolate these low energy 
interaction properties to higher energies. Consequently they provide different expectations 
for (X max ) and <r(A max ). The first aim of this paper is to discuss how the mean value of the 
depth of shower maximum and its dispersion can be used to interpret mass composition even 
in the presence of uncertainties in the hadronic interaction modeling. 

Furthermore, we discuss the different roles of the two observables, (A max ) and cr(A max ), 
with respect to mass composition. In the interpretation of data they are often used as 
different, and independent, aspects of the same phenomenon. However it is not true to say 
that both parameters reflect the cosmic ray composition to the same extent. According to 
the superposition model [2] (A max ) is linear in (In A) and therefore it actually measures mass 
composition for both pure and mixed compositions. But, we will show that the behaviour of 
cr(X max ) is more complex to interpret as there is no one-to-one correspondence between its 
value and a given mean log mass. Only in the case of pure composition is this correspondence 
unique. 

In this paper we refine the analysis method originally proposed by Linsley [3, 4] and 
apply it to the Auger data. The Pierre Auger Collaboration has published results on the 
mean and dispersion of the A max distribution at energies above 10 18 eV [5, 6]. In this work 
we apply the proposed method to convert those observables to the first moments of the log 
mass distribution, namely (In .A) and crf nA - 

The paper is organized as follows. In Section 2 we discuss the parameterization for 
(A max ) and 0"(A max ). In Sec. 3 we test the method with shower simulations assuming 
different mass distributions. Sec. 4 describes the application of the method to data. The 
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discussion of the results and the conclusions follow in sections 5 and 6 respectively. The 
details of the parameterization and the best fit values for the hadronic interaction models 
are summarized in Appendix A. 



2 A method to interpret (A" max ) and a(X mayL ) 

The interpretation of (A max ) and 0"(A max ) can be simplified by making use of an analysis 
method based on the generalized Heitler model of extensive air showers [7]. In this context 
(A max ) is a linear function of the logarithm of the shower energy per nucleon: 

(X^)=X + D]og 10 (J^\ , (2.1) 

where Ao is the mean depth of proton showers at energy Eq and D is the elongation 
rate [8-10], i.e., the change of (A max ) per decade of energy. The High Energy hadronic 
interaction models used in this work are EPOS 1.99 [11], Sibyll 2.1 [12], QGSJet 01 [13] and 
QGSJet II [14]. Simulated data show that eq. (2.1) gives a fair description of EPOS and 
Sibyll results in the full range of interest for this work, 10 18 to 10 20 eV, but does not reproduce 
accurately QGSJet models. For this reason we generalize the original representation as: 

(A max ) = A + D log 10 (^) +£\nA + 5 In A log 10 (J^j , (2.2) 

where the parameters £ and 6 are expected to be zero if the model predictions are compatible 
with the superposition result (2.1). 

For nuclei of the same mass A one expects the shower maximum to be on average: 

(A max ) = (A max ) p + f E In A , (2.3) 

and its dispersion to be only influenced by shower-to-shower fluctuations: 

<r 2 (A max ) = <7 s 2 h (lnA) . (2.4) 

Here (A max ) p denotes the mean depth at maximum of proton showers, as obtained from either 
eq. (2.1) or (2.2), and o" 2 h (lnA) is the A max variance for mass A, c 2 h (ln^4) = o" 2 (A max | In A). 
The energy dependent parameter appearing in (2.3) is: 

'* = «-5no +,k *»(f)- <2 ' 5) 

The values of the parameters Ao, D, £, 5 depend on the specific hadronic interaction model. In 
this work they are obtained from CONEX [15] shower simulations as described in Appendix A. 

In the case of a mixed composition at the top of the atmosphere, the mean and variance 
of A max depend on the lnyl distribution. There are two independent sources of fluctuations: 
the intrinsic shower-to-shower fluctuations and the InA dispersion arising from the mass 
distribution. The first term gives rise to (cr 2 h ), the average variance of A max weighted 

( d{X max ) \ 2 2 



according to the InA distribution. The second contribution can be written as y d\nA ) °iiM 

where crf nA is the variance of the In A distribution. We can finally write for the two profile 
observables: 

(A max ) = (X max ) p + f E (In A) (2.6) 
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0- 2 (^max) = (cr s 2 h ) + /| af nA . (2.7) 

The two equations depend on energy through the parameters but also via (c 2 h ) and the 
possible dependence of the two moments of the In A distribution. 

To obtain an explicit expression for (<7 2 h ) we need a parameterization for a 2 h (lnA). We 
assume a quadratic law in lnj4: 

a s 2 h (lnA) = o 2 p [l + a\nA + b{\nA) 2 } , (2.8) 

where a 2 is the X miiX variance for proton showers. The evolution of a 2 h (lnA) with energy is 
included in a 2 and the parameter a: 



°l = Po + Pi logio ( TT ) + P2 



E 



E 



log 



10 



E 



E: 



2 



The parameters po, p±, p2, ao, a±, b depend on hadronic interactions: the values used in the 
paper are given in Appendix A. 

Using measurements of (A max ) and <r(A" max ), equations (2.6) and (2.7) can be inverted 
to get the first two moments of the In A distribution. From eq. (2.6) one gets: 

(lnA) = (*™*)-(*™*)p , (2.io) 

JE 

Averaging eq. (2.8) on ln^4 one obtains: 

(a 2 h ) = a 2 p [l + a(lnA)+b((lnA) 2 )] . (2.11) 
Substituting in eq. (2.7) we get: 

a 2 (X mSLX ) = a 2 [l + a(hxA) + b((lnA) 2 }] + f 2 a 2 nA . (2.12) 
But by definition ((In A) 2 ) = af nA + (In A) 2 . Solving in af nA one finally obtains: 

- ^TTI ' (2 ' 13) 

Equations (2.10) and (2.13) are the key tools used throughout this work for interpreting 
Pierre Auger Observatory data in terms of mass composition and assessing the validity of 
available hadronic interaction models. 



3 Testing the method with simulation 

Equations (2.6) and (2.7) can be tested with simulations. They contain parameters depending 
on the hadronic interaction properties and on the mass distribution of nuclei. The mass 
distribution of nuclei refers to those nuclei hitting the Earth's atmosphere: it does not matter 
what is the source of the mass dispersion, either a mixed composition at injection or the 
dispersion caused by propagation. So, in order to test the method we will simply use different 
test distributions of the masses at the top of the atmosphere. 

For this purpose we have chosen three different mass distributions: 

1. A distribution uniform in ln^4 from ln(l) to ln(56) and independent of energy. The 
values of {In A) and a\ n A are respectively 2.01 and 1.16. 
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2. A Gaussian In A distribution with (In A) increasing linearly with logE from In (4) at 
10 18 eV to ln(14) at 10 20 eV and <r\ n A = 0.75 independent of energy. The Gaussian is 
truncated to less than 2 sigmas to avoid unphysical mass values. In this case the InA 
dispersion is fixed and equal to 0.66 but (In .A) varies with energy. 

3. Two masses, H and Fe, with proton fraction H/(H + Fe) decreasing linearly with logE 
from 1 at 10 18 eV to at 10 20 eV. In this case, both (In A) and <J\nA vary with energy. 






Figure 1. (X max ) and <r(X max ) as a function of log 10 (-E/eV) for three different mass distribution 
hypotheses (see text). Full circles are calculated from the resulting X max distributions from the 
CONEX simulations. Sibyll 2.1 has been chosen for hadronic interactions. The dashed lines show 
equations (2.6) for (X max ) and (2.7) for er(X max ). The dot-dashed line refers to the contribution of 
the first term in (2.7). 

Figure 1 shows the result of the test for the three mass distribution hypotheses. To 
generate the X max distributions we have used CONEX [15] showers with Sibyll 2.1 [12] as 
the hadronic interaction model. These distributions do not include detector effects. For each 
test mass hypothesis, the mean and RMS are retrieved from the resulting X max distribution 
obtained from the simulations. These are shown as full circles, (X max ) and cr(X max ) in left 
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and right panels respectively. The dashed lines are calculated using equations (2.6) and (2.7) 
for the three different mass hypotheses by using only the first two moments (In .A) and a\ Q A- 
One can see that, despite the simple assumptions made, good agreement is achieved for 
all the three mass distributions. The dot-dashed line refers to the contribution of the first 
term in eq. (2.7). The comparison between the two lines (dashed vs. dot-dashed) highlights 
how different the interpretation of cr(X max ) data can be if one does not take into account the 
mass dispersion term. 

The inverse equations (2.10) and (2.13) have also been tested using Monte Carlo sim- 
ulation. In this case (In A) and crf nA have been obtained as a function of log 10 (-E/eV) 
directly from the input mass distributions. These values are shown as full circles in Figure 
2. The (X max ) and cr(X max ) retrieved from the corresponding X max distributions are used 
in equations (2.10) and (2.13) to get (In A) and cr^ A . These are shown in Fig. 2 as dashed 
lines. Also in this case, the comparison is quite successful. 
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Figure 2. (In A) and crf nA as a function of log 10 (£'/eV) for three different mass distribution 
hypotheses. Sibyll 2.1 is the hadronic interaction model. Full circles refer to the values obtained 
directly from the input mass distributions. The dashed lines show (In A) and af nA calculated using 
equations (2.10) and (2.13). The dotted lines refer to the calculation of the same variables using the 
parameterization for QGSJct II in (2.10) and (2.13). 
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The simulated data sample can also be used to estimate the systematic uncertainty 
in the calculation of the moments of the X max (In A) distribution induced by the missing 
knowledge of the hadronic interaction mechanism. This study is pursued using simulated 
showers generated with a given model together with parameters of another model in equations 
(2.6) and (2.7) for the profile variables, and (2.10) and (2.13) for the log mass variables. An 
example of this procedure is shown in Fig. 2 where the dotted lines show the calculation 
with the parameters of QGSJet II and the full circles refer to data simulated with Sibyll 
2.1. As a summary of these cross- model checks, we find mean absolute deviations of 4 to 
27 g cm" 2 for (A" max ) and 1 to 5.4 g cm -2 for cr(X max ), where the maximum deviations are 
obtained crossing EPOS with QGSjetll. The same study done for the moments of the log 
mass distribution gives mean absolute deviations of 0.2 to 1.2 for (In A) and 0.02 to 0.5 for 
crf nA . In this case the maximum values refer to EPOS vs. QGSJet 01 for the first moment 
and QGSJet II vs. QGSJet 01 for the second. 



4 Application to data 

At ultra-high energies, shower development can be directly measured using fluorescence and 
Cherenkov light profiles. Mean X m3X data as a function of energy are available from Fly's 
Eye [16], HiRes [17, 18], Auger [5], Yakutsk [19] and Telescope Array [20]. (A~ max ) data were 
complemented with fluctuation measurements as early as 1980s (see e.g. [21] and references 
therein) but only recently have precise optical detector measurements become available [5, 
18, 19]. 

The Pierre Auger Collaboration has published results on the mean and dispersion of 
the X max distribution at energies above 10 18 eV [5]. Here we apply the method presented in 
this work to an updated dataset available in [6, 22]. These data are shown in Figure 3. 
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Figure 3. (X max ) (left) and <r(X max ) (right) as a function of log 10 (-B/eV) from Pierre Auger 
Observatory data [6, 22]. Data (full circles) are shown with statistical errors. Systematic uncertainties 
are represented as bands. 



In the Auger analysis [5] , the events are selected using fiducial volume cuts based on the 
shower geometry. This ensures that the viewable X m3iX range for each shower is large enough 
to accommodate the full X max distribution. Also, the detector resolution is accounted for 
by subtracting in quadrature its contribution to the measured dispersion. This allows the 
direct conversion to the moments of the In A distribution using equations (2.10) and (2.13) 
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without the need of more complex treatment, such as is required in the presence of acceptance 
biases [23, 24]. 

The moments of the log mass distribution, (In A) and crf nA , as obtained using equations 
(2.10) and (2.13), are shown (full circles) as a function of log± Q (E/eV) in Figures 4 and 
5 respectively. Error bars show the statistical errors obtained from the propagation of 
data errors and the errors of the fitted parameters. Shaded bands are the systematic 
uncertainties obtained by summing in quadrature the different individual contributions. 
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Figure 4. {in A) as a function of log 10 (£'/eV) obtained from Auger data [22] are shown as full circles 
for different hadronic interaction models. Error bars show statistical errors. The shaded areas refer to 
systematic uncertainties obtained by summing in quadrature the systematic uncertainties on (X max ) 
and cr(^fmax) data points and on the FD energy scale. 

The systematic uncertainties on (X max ) and cr(X max ) data points have different sources: 
calibration, atmospheric conditions, reconstruction and event selection [5]. Another source 
of systematics is related to the uncertainty of the FD energy scale [25], 22 %, which induces 
an uncertainty in (In A) and crf nA via the parameters of the models. All these uncertainties 
contribute approximately at the same level and independently of energy. The figures show 
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Figure 5. (jf nA as a function of log 10 (i?/eV) obtained from Auger data [22] are shown as full circles 
for different hadronic interaction models. Error bars show statistical errors. The shaded area refers 
to systematic uncertainties as in Fig. 4. The lower limit of allowed &f nA is shown by the exclusion 
line. The upper limit (4.05) is just above the maximum of the vertical axis. 



the results for the moments of the log mass distribution for EPOS 1.99 [11], Sibyll 2.1 [12], 
QGSJet 01 [13] and QGSJet II [14]. 

Despite the uncertainties and the different mass offsets of the models, the overall features 
are similar in all the cases. So far as the energy dependence is concerned, the data imply 
an increasing (ln^4) above 10 18 ' 3 eV from light to intermediate masses and a decreasing cr^ nA 
over the whole energy range. 

Looking more specifically to the different hadronic models we notice a slight change 
in the log mass scale. The highest masses are obtained for EPOS 1.99. Sibyll 2.1 and 
QGSJet II show intermediate values, whereas the lowest masses are obtained for QGSJet 01. 
In particular at log 10 (£'/eV) = 18.25 the mean log mass, (In A), is 1.10, 0.70, 0.60 and 0.12 
respectively for EPOS 1.99, Sibyll 2.1, QGSJet II and QGSJet 01 with statistical errors of 
about 0.08 and systematic uncertainty of about 0.6. The Pierre Auger Collaboration has 
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recently published the measurement of the proton-air cross section for the energy interval 
10 18 to 10 18,5 eV [26]. That measurement is done using the showers with A max > 768 g cm~ 2 , 
corresponding to 20% of the total X max distribution. Even in the most unfavourable case, 
(the (In A) and of nA predicted by EPOS), one finds that several realizations obtained from 
the allowed (In A) and crf nA have enough protons in the most deeply penetrating showers to 
fulfill the selection criteria adopted in the Auger analysis. 

Whereas (In A) always has valid values (apart a small region which crosses (In A) = for 
QGSJet 01), there are wide energy intervals where crf nA is negative. Considering eq. (2.13) 
one can see that these values occur for energies where the shower fluctuations corresponding 
to the mean log mass exceed the measured A max fluctuations. Figure 5 shows that crf nA data 
points are within the allowed physical region only for EPOS 1.99 and Sibyll 2.1. They are 
partly outside for QGSJet II, and completely outside for QGSJet 01. However the current 
systematic uncertainties do not allow one to establish stringent tests to the models. 

The method presented in this work shows that the Pierre Auger Observatory data can 
confront hadronic physics models provided that future developments in the shower data 
analysis reduce systematics. By shrinking the shaded bands in Figure 5 it will be possible to 
constrain those models. 

5 Discussion 

The importance of the combined study of the mean values and fluctuations of mass dependent 
observables has been addressed by several authors [3, 4, 21, 27, 28]. In particular, Linsley [4] 
showed that a combined analysis of the mean and the variance of In A can provide a useful 
representation of the mass transition (if any) to be found in shower profile data. In fact, pos- 
sible transitions are constrained to a limited region of the ((In A), of nA ) plane. More recently 
a similar study using the (A max )-o"(A max ) correlation 1 reached a similar conclusion [29]. 

Converting X max data to In A variables, as described in Sec. 2, one can plot Pierre 
Auger Observatory data in the ((In A), of nA ) plane. Since this procedure depends on the 
hadronic model, one gets a plot for each model as shown in Figure 6. Data points are shown 
as full circles with size increasing in proportion to logE. The error bars are tilted because of 
correlations arising from equations (2.10) and (2.13) and represent the principal axes of the 
statistical error ellipses. The solid lines show the systematic uncertainties. The same figure 
shows the region allowed for mass compositions. The contour of this region (gray thick line) 
is generated by mixing neighbouring nuclei in the lower edge and extreme nuclei (protons 
and iron) in the upper edge. Each of these mixings is an arch shaped line in the ((In A), 

a inA) P lane - 

Figure 6 shows that the Auger data lie outside the allowed boundaries for part of the 
energy range in some of the models. As noted previously, systematic uncertainties are still 
large and thus prevent us from more definite conclusions. However the energy evolution is 
common to all models suggesting that the average mass increases with decreasing log mass 
dispersion. This behaviour might imply astrophysical consequences. 

In fact there are only a few possibilities for extragalactic source models to produce 
compositions with small log mass dispersion at the Earth. Protons can traverse their path 
from sources to the Earth without mass dispersion, but this case is excluded by Pierre Auger 
Observatory data at the highest energies. 

1 In this case the dependence on hadronic models has been accounted for by subtracting the corresponding 
observables predicted by the models for iron. 
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Figure 6. Pierre Auger data in the ((In A), <J? nA ) plane for different hadronic interaction models. 
Data points are shown as full circles with statistical errors. The marker sizes increase with the 
logarithm of the energy. Systematic uncertainties are shown as solid lines. The gray thick line shows 
the contour of the (In A) and <J\ nA values allowed for nuclear compositions. 



Nuclei originating from nearby sources (j$ 100 Mpc) might be detected with small mass 
dispersion. For these sources, propagation does not degrade mass and energy so the spectrum 
and composition reflect closely their values at injection. But, if sources are distributed 
uniformly, distant sources induce natural mass dispersions. Small In A dispersions are possible 
only when there is small observed mass mixing so that, at each energy, only nuclei with a 
small spread in masses are present. This corresponds to the low-uj^ edge of the contour of 
the allowed region in the ((In A), <rf nA ) plane. 

Protons originating by the photo-disintegration of nuclei are the main source of mass 
dispersion because they populate each energy region. The possible end of the injection 
spectrum based on a rigidity-dependent mechanism can reduce the proton component at 
high energies, thus producing a reduction of the mass dispersion at the highest energies. A 
complete study of source models under several hypotheses is required to study all the source 
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parameters that limit the mass dispersion in the propagation of extragalactic cosmic rays. 
Recent studies, see e.g. [30, 31], based on the assumption of a uniform source distribution, 
have shown that the Auger composition results, when combined with the energy spectrum, 
require hard injection spectra (i.e. index < 2) with low energy cutoffs and the possible 
presence of local sources. 

6 Conclusions 

In this work we presented a method for interpreting (A max ) and cr(X max ) in terms of mass 
composition. The method is based on an extension of the Heitler model of extensive air 
showers. The parameterization given in equations (2.6) and (2.7) expresses those two profile 
observables as a linear combination of the first two moments of the log mass distribution, 
(In A) and cr^ nA , and of the mean shower fluctuations. 

We first note that the method provides an effective key to the interpretation of data. The 
energy dependences of (A" max ) and cr(X max ) are sometimes considered as different expressions 
of the same physical features, e.g. an increase or decrease of the mean log mass. However 
their different meanings can be easily understood by looking at the dependence on the mass 
variables. At a fixed energy (A" max ) is only function of (In A); therefore, it only carries 
information of the average composition. However, cr(X max ) cannot be interpreted as a 
measure of the average composition since it is also affected by the log mass dispersion. 
Similarly, the inference of hadronic interaction properties from cr(X max ) can be wrong unless 
the mass dispersion term (oc crf nA ) is negligible. The parameter a\ n A represents the dispersion 
of the masses as they hit the Earth atmosphere. It reflects not only the spread of nuclear 
masses at the sources but also the modifications that occur during their propagation to the 
Earth. 

The method has been succesfully tested, with the simulation of different mass distri- 
butions in the energy interval from 10 18 to 10 20 eV showing the robustness of the param- 
eterization. We have applied the method to the Pierre Auger Observatory A max data to 
get the first two moments of the In A distribution. The outcome relies on the choice of a 
hadronic interaction model to set the parameters and the appropriate shower fluctuations. 
Four models have been used, EPOS 1.99, Sibyll 2.1, QGSJet 01 and QGSJet II, and the 
corresponding moments of the log mass distribution have been obtained as a function of 
energy. Despite the differences in the chosen models, the overall features are quite similar. 
In particular we find an increasing (In A) above 10 18,3 eV from light to intermediate masses 
and a decreasing crf nA over the whole energy range, while the mean log mass scale changes 
with hadronic models. 

The results presented in this paper show the capability of the method to infer important 
features of the mass distribution of the UHECR nuclei. This is a remarkable outcome with 
respect to the study of source scenarios and propagation. In fact we do not only access the 
average mass, but also the mass dispersion. While a pure proton beam at the sources is not 
changed by propagation, nuclei should increase the mass dispersion in their path towards the 
Earth. The Auger results seem to imply either close-by sources or hard spectral indices, if 
the energy evolution of the present hadronic interaction models can be trusted. 

The proposed method can also be used as a tool to investigate the validity of hadronic 
interaction models. In particular it has been shown that the intrinsic shower fluctuations 
are sometimes larger than the measured A max dispersions. This happens in different en- 
ergy intervals for the different models. At the highest energies, all models approach the 
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lower boundary, and some of them enter the unphysical region, but the current systematic 
uncertainties prevent us from confidently rejecting any model. Provided that systematic 
uncertainties can be reduced in future data analysis, the method can be used to constrain 
hadronic interaction models. The addition of new measurements, such as the muon content 
of EAS [32, 33], may allow us to place stronger bounds to the models. 

A Parameterization of the shower mean depth and its fluctuations 

The shower code chosen for this work is CONEX [15]. CONEX is a hybrid simulation code 
that is suited for fast one-dimensional simulations of shower profiles, including fluctuations. 
It combines Monte Carlo simulation of high energy interactions with a fast numerical solution 
of cascade equations for the resulting distributions of secondary particles. In our CONEX 
simulation we used the default energy thresholds settings of version v2r3.1 2 . 

The parameters Xq, D\, £ and 5 used in equations (2.6) and (2.7) have been obtained 
by fitting CONEX showers for four different primaries (H, He, N and Fe) in nine energy bins 
of width A log 10 (-E/eV) = 0.25 ranging from 10 18 to 10 20 eV, and for all the hadronic models 
used in this work: EPOS 1.99 [11], Sibyll 2.1 [12], QGSJet 01 [13] and QGSJet II [14]. 
In total, about 25,000 showers have been used for each energy bin and for each hadronic 
model. The fit procedure always converges with mean (maximum) (A" max ) residuals from the 
simulated data of about 1 (3) g cm -2 for all the models. The best fit values are reported in 
Table 1 with their errors. 



parameter 


EPOS 1.99 


Sibyll 2.1 


QGSJet 01 


QGSJet II 


X 


809.7 ± 0.3 


795.1 ± 0.3 


774.2 ± 0.3 


781.8 ± 0.3 


D 


62.2 ± 0.5 


57.7 ± 0.5 


49.7 ± 0.5 


45.8 ± 0.5 




0.78 ± 0.24 


-0.04 ± 0.24 


-0.30 ± 0.24 


-1.13 ± 0.24 


<5 


0.08 ± 0.21 


-0.04 ± 0.21 


1.92 ± 0.21 


1.71 ± 0.21 



Table 1. Parameters of formulae (2.6) and (2.7) for different hadronic interaction models setting Eq 
= 10 19 cV. The values are obtained fitting the mean X max for showers generated for four different 
primaries H, He, N and Fe. Statistical error obtained from the fit are also given. All values are 
expressed in g cm~ 2 . 

Shower variances have been fitted using the parameterization given in equations (2.8) 
and (2.9) and the same simulated data set described above. The mean (maximum) <r(A max ) 
residuals from the simulated data are about 1 (3) g cm~ 2 for all the models. The best fit 
parameters are given in Table 2 with their errors. 
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parameter 


EPOS 1.99 


Sibyll 2.1 


QGSJet 01 


QGSJet II 


po x (g~ 2 cm 4 ) 


3279 ± 51 


2785 ± 46 


3852 ± 55 


3163 ± 49 


Pi x (g~ 2 cm 4 ) 


-47 ± 66 


-364 ± 58 


-274 ± 70 


-237 ± 61 


P2 x (g cm J 


228 ± 108 


152 ± 93 


169 ± 116 


60 ± 100 


a 


-0.461 ± 0.006 


-0.368 ± 0.008 


-0.451 ± 0.006 


-0.386 ± 0.007 


a\ 


-0.0041 ± 0.0016 


-0.0049 ± 0.0023 


-0.0020 ± 0.0016 


-0.0006 ± 0.0021 


b 


0.059 ± 0.002 


0.039 ± 0.002 


0.057 ± 0.001 


0.043 ± 0.002 



Table 2. Parameters of formulae (2.8) and (2.9) for different hadronic interaction models setting E$ 
= 10 19 eV. The values are obtained fitting cr 2 (X max ) for showers generated for four different primaries 
H, He, N and Fe. The statistical errors obtained from the fit are also given. 
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